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i ' Abstract 

A new package to simulate the formation of electrical pulses in segmented true- 
coaxial high purity germanium detectors is presented. The computation of the elec- 
tric field and weighting potentials inside the detector as well as of the trajectories of 
the charge carriers is described. In addition, the treatment of bandwidth limitations 
and noise are discussed. Comparison of simulated to measured pulses, obtained from 
an 18-fold segmented detector operated inside a cryogenic test facility, are presented. 
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1 Introduction 
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High purity germanium detectors, HPGeDs, are used in a wide variety of 
applications in particle and nuclear physics [T|l2] . The analysis of their pulse 
shapes, PSA, can be a powerful tool to determine event topologies, often a key 
to distinguish between signal and background events [3 j4T5l6| . In particular, 
the separation of events with energy depositions in one or multiple locations 
inside an HPGeD, single-site or multiple-site events, is of interest. Previous 
studies were often limited by the availability of clean event samples to train 
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neural nets or to evaluate PSA [5f6] . Therefore, simulated events reflecting a 
detailed understanding of the underlying processes are extremely valuable. 

The interactions of radiation (a, fl, 7, n,p, etc.) inside a semiconductor create 
electron-hole pairs. These charge carriers are separated by the electric field 
inside an HPGeD and drift towards the electrodes inducing time dependent 
charges. In most applications, charge sensitive amplifiers are used and the 
resulting pulses are sampled and digitized at a given frequency. 

A new pulse shape simulation, PSS, package was developed. The geometry 
of an 18-fold, (60, 3z), segmented true-coaxial HPGeD [7] is implemented as 
default. The inner mantle has a radius of 5 mm, the outer of 37.5 mm. The 
length is 70 mm. This allows the comparison to data obtained from such a 
detector which was developed in connection to the GERDA [8] experiment. 
However, the simulation package itself is general and the geometry can easily 
be modified to accommodate any kind of germanium detector. 



2 Simulation procedure 

The package P|TU] is separated into two parts, the calculation of the static 
properties of the HPGeD followed by the event by event simulation of the 
pulse development. The tabulation of the static properties is done once at the 
beginning of the simulation. It includes the electric field and the weighting 
potentials inside the detector. 

The event by event simulation comprises the event topology and the develop- 
ment of the pulses in time-steps: 

(1) Simulation of the interactions of particles with germanium. The output 
are individual positions with energy depositions, hits; 

(2) Clustering of hits. According to the requirements of the simulation, hits 
can be clustered 0. The position of the cluster is the bary-center of the 
original hits. The energy of the cluster, E^, is the sum of the energies of 
the original hits. 

(3) Simulation of the drift of charge carriers in time-steps. For each cluster, 
only one point-like charge is considered. At each time-step the velocity 
of the charge is calculated using the electric field as tabulated; 

(4) Calculation of the charges induced in the electrodes after each step; 

(5) Simulation of experimental effects such as noise, bandwidth limitation, 
shaping times, etc.. 

1 Hits closer to each other than the best hypothetical radial resolution of the de- 
tector, 1 mm, calculated as the speed of the charge carriers multiplied by the time 
resolution of the electronics, are clustered by default. 
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The first step is done using the GEANT4 [TTp2] based simulation package 
MaGe[13], jointly developed by the GERDA and Majorana collaborations. 
The good agreement between MaGe and data obtained with segmented de- 
tectors was demonstrated previously [3]. 

All calculations are performed in cylindrical coordinates, r, and z, with the 
origin at the center of the detector. 



3 Electric field and weighting potentials 

The electric field, E(x), at a position x = (r,<f>,z) inside an HPGeD depends 
on the geometry of the detector, the bias voltage, Vu a s applied and the density 
of electrically active impurities, pi mp (x). It is considered static, calculated only 
once and tabulated for a reasonably spaced grid 0. The electric field at any 
position is calculated by interpolation. 

For simple cases like a constant Pi mp {x) , E(x) can be calculated analytically 
solving Poisson's equation \7 -E(x) = where p(x) is the space charge den- 
sity e is the dielectric constant. However, the electric field is always calculated 
numerically using the potential ip(x) and the boundary conditions. 

The boundary conditions for <p depend on V^as- For the true-coaxial n-type 
HPGeDs considered here, the potential is fixed to Vu aa on the inner mantle, 
zero on the outer mantle and it floats on the end surfaces. The density Pi mp (x) 
determines p(x) and thereby affects the electric field inside an HPGeD. In the 
case of constant pi mp , the electric field only depends on r. For this case, the 
output of the numerical calculation based on successive over-relaxation was 
tested against the numerical solution pU]. The deviations for realistic values 
of Pim P (x) were found to be less than 0.6% throughout the detector. 

Figure [1] shows the calculated strength of the electric field as a function of r for 
Vbias = 3kV and varying constant pi mp . This range of pi mp values is realistic 
for HPGeDs. For lower pi mp , the field is strong throughout the detector. For 
higher pi mp , the field is insufficient at small radii; a higher Vu as is needed. 

As Pim P (x) can vary up to a factor three between the two ends of a detector and 
might depend on r, it is possible to specify a three dimensional distribution 
for p imp {x). 

The influence of any electrode on a given space-point inside the HPGeD is 
characterized by its weighting potential. This is defined according to Shockley- 
Ramo's theorem [23l24}25j as the solution of Poisson's equation for the bound- 

2 The default is a grid with a 1 mm spacing in r and z and a 2 degree spacing in <fi. 



3 




E 

:> 

LU250- 



P 
P 
P 



0.3- 10 10 cm- 3 
0.5- 10 10 cm- 3 
0.7- 10 10 cm- 3 
0.9- 10 10 cm- 3 
1.1-10 10 cm- 3 
1.3- 10 10 cm- 3 



200- 



100-V 



150 — 




— v 



50^-. 




i u» 

10 



J L 




15 



20 



30 



Radius [mm] 



Fig. 1 . Strength of the electric field as a function of the cylindrical coordinate r for 
constant impurity densities between 0.3 and 1.3 x 10 _10 /cm 3 . 

ary conditions that the potential on the electrode of interest equals unity and 
the potentials on all other electrodes equal zero. All weighting potentials are 
treated numerically like the electric field. 



4 Drift of charge carriers 

The drift of the charge carriers is calculated in time-steps, At, for each point 
representing a cluster of hits as defined in Section [2j At each step the velocity 
is calculated using the electric field as tabulated. Between grid points a sim- 
ple linear extrapolation is used for the electric field. Two different numerical 
methods are implemented to calculate the trajectories, the Euler method and 
the 4 th Runge-Kutta method. The former is less computer time intensive, but 
is also less precise. However, for step-sizes At < 1 ns, the resulting trajectories 
do not differ significantly [9] . 

The drift velocity of the charge carriers, v e (x) for electrons and Vh(x) for holes, 
is calculated for each step using the electric field, E( tabulated: 



where fi e / h (x) is the mobility pQ of electrons and holes, respectively. 
Whether fi e /h(x) depend on the relative position of x to the crystal axes is 



Ve/h(x) = II e/h (x)E(x), 



(1) 
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determined by the temperatures of the germanium crystal, T crysta i, and of the 
charge carriers, T cc . If T crysta i rs T cc , fi e / h (x) become numbers, and the 
drift velocities are directly proportional to the electric field. However, HPGeDs 
are normally operated at temperatures around 100 K and T crysta i « T cc . In 
this case, the mobility becomes a complex tensor. Therefore, the simulation 
has to take the crystal structure into account and the calculation of v e /h{x) 
becomes quite involved. The drift of the charge carriers is not any longer 
parallel to the electric field [H] for all x. 

The axes of the face-centered cubic structure of germanium crystals are char- 
acterized by Miller indices jTS] . Cylindrical germanium detectors are produced 
with their z axis aligned to the crystal axis (001) [16] . This is implemented in 
the simulation. The axes (100) and (110) are in the r0-plane. Their relative 
position to the segment boundaries can be chosen freely, i.e. adjusted to the 
detector to be simulated. 

The tensors fi e /h{x) are unfortunately not known. Special measurements along 
the crystallographic axes are used and parameterized. The results are com- 
bined to get results for any point in the crystal. 

If the electric field is parallel to any of the principal crystallographic axes, the 
charge carriers will drift along this axis, because of the symmetric structure 
of the germanium crystal [H][TT] . The measured drift velocities [T8fl9f20] . 
v e /hi along the axes (100) and (111) with E(x) || (100) for x on (100) and 
E(x) || (111) for x on (111), respectively, are parametrized well [TfT9] for 
E = \E(x)\ < 300 V/mm by 

Ve/h [1 + (M£l)/3]l//3' 1 ' 

where fj^/ h , E Q and (3 are parameters determined by fitting. The parameters 

H Q e , h represent the linear relation between v and E at large T crysta i and low 
E. The parameters E Q and /3 are used to model the deviation from this linear 
relation at low lattice temperature and high electric fields. 

The set of parameters given in |18j and validated for the drift of electrons [2TJ 
was used in the simulation presented here. The values are listed in Table [TJ 

The drift velocities in any direction can be derived from the velocities along 
the (111) and (100) axes. The model used for the electron drift [19] is based 
on the idea that the conduction band in a germanium crystal reaches its min- 
imal potential in regions around the four equivalent (111) axes. Free electrons 
effectively only populate these regions. The probability density of free elec- 
trons in other regions can be ignored. The resulting drift velocities along the 



5 



Table 1 

Parameters for the experimental drift velocities in the (111) and (100) directions 
used in the simulation. 
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Fig. 2. Drift velocities of (a) electrons and (b) holes along the principal crystal 
axes as function of the electric field. Velocities along the axes (111) and (100) are 
calculated according to Eq. [2j For the calculation along (111), see text. 

(110) axis are shown together with the velocities for the (111) and (100) axes 
in Fig. E{a). 

The model used to calculate the hole drift velocities [IB] in any direction is 
based on the idea that only the "heavy hole valence band" [22] is responsible 
for the anisotropy of the mobility. All other effects are neglected. A parameter- 
ization [18] is used to calculate the three components of the hole drift velocity 
v at any position. The results for the (110) axis together with the values for 
the (111) and (100) axes are shown in Fig. |2{b). 



5 Trajectories 

The results presented here were obtained with the Runge-Kutta method. 
Figure [3] shows the trajectories of electrons and holes created at the outer 
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Fig. 3. Trajectories of (left) electrons and (right) holes in the (x, y) = (r, <p) plane 
for constant pi mp . Electrons drift inwards and holes outwards. 

and inner mantle, respectively. The step size was set to 1 ns and the total 
time simulated was 400 ns, Vu a s = 3000 V and p imp = 0.62 • 10 10 was constant 
through the volume. The trajectories are bent due to the crystal structure; 
this is called transverse anisotropy. The charge carriers are slower in the (110) 
than in the (100) direction; this is called longitudinal anisotropy. As holes 
are slower than electrons, the longitudinal anisotropy can be seen very clearly 
along (110) where the holes do not reach the outer mantle within 400 ns. 

The transverse anisotropy is also demonstrated for electrons in fig. |H It shows 
the displacement A0 at the end of the drift and the total time needed for the 
drift as functions of the azimuthal angle at the start. Along the crystal axes 
(110) and (100), A0 = 0. The drift time is minimal at (100) and maximal at 
(110). 



6 Charges induced on electrodes 



The electric signals are induced in the electrodes of a detector by the cu- 
mulative influence of moving electrons and holes. Shockley-Ramo's Theorem 
[2"3"|24|I25] is used to calculate the time development of the induced charge Q{t) 
or current I(t) in each electrode: 

Q(t) = -Q ■ [<fw{Xh(t)) - <Pw(%e{t))], (3) 

I(t) = Q ■ [E w (x h (t)) ■ v h (t) - E w (x e (t)) ■ v e {t)}, (4) 

where Q is the electric charge carried by electrons or holes, x e /h(t) and v e /h(t) 
are the position and velocity vectors of electrons/holes as a function of time, 
and <f w and E w are the weighting potentials and weighting fields. 
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Fig. 4. (a) Displacement Acfi between starting point on the outer mantle and the 
end point on the inner mantle of the detector as a function of the azimuthal angle 
4> at the start. (b) Total drift-time as a function of 4>. 



Figure^ shows the weighting potential of a segment, B, with the indication 
of a photon interaction. The drift trajectories of holes and electrons are also 
indicated. Figures 0(b) and (c) show the charges and currents induced on the 
electrodes of the segment with the interaction and its neighbors. The pulses 
induced in the neighboring segments are called mirror pulses. The amplitude 
of the mirror pulse induced in segment A is larger than the one in segment C 
since the trajectory is closer to segment A than C. 

The charge pulses in segment B and the core show clear kinks at the time when 
the inwards drifting electrons reach the core electrode, i.e. the inner mantle. 
At that time the current pulses show maxima. The charge mirror-pulses do not 
show clear features at this time; the current mirror-pulses show more distinct 
features. 

The position in r of an energy deposit is clearly reflected in the resulting pulses 
in the segment and core electrodes. This is demonstrated in Figs. |6] and CD 
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Fig. 5. (a) weighting potential of segment B together with an indication of a 7 
interaction, (b) Simulated charge and (c) current pulses induced in segments A, B 
and C and the core. 

7 Influence of electronics 

The most widely used amplification scheme is based on charge sensitive de- 
vices. The bandwidth is limited by the amplifiers and by the transmission 
cables. In order to account for the limited bandwidth and the decay time of 
the amplifier, the simulated pulses can be convoluted with a transfer function 
representing the readout electronics. In addition, Gaussian noise can be added 
to represent the noise of the system. 

Figure [8] demonstrates the modification of a simulated pulse after folding in 
a decay time of 5 /is and a cut-off in bandwidth was 10 MHz. The noise level 
was set to 5 % of the pulse amplitude. These noise and band-width values are 
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Fig. 6. Simulated charge pulses in the segment electrode for charges deposited at 
different radius, r. 
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Fig. 7. Simulated charge pulses in the core electrode for charges deposited at different 
radius, r. 

chosen to clearly demonstrates the effects. 

The decay-time of the preamplifier was chosen significantly larger than the 
length of the pulse. Therefore, mainly the plateau is affected. The cut-off in 
bandwidth smoothes the features of the pulse. The kink clearly visible at 
t = 150 ns in Figure IU(b) is washed out significantly. The noise applied for 
Figure ffl^d) makes the kink almost invisible. 
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Fig. 8. Pulses at different stages of the simulation: (a)Raw pulse, (b) after folding 
in in the decay time of the signal, (c) after convolution with the bandwidth and (d) 
after adding noise. 



8 Uncertainties of the simulation 



The output of the simulation depends critically on the input parameters used. 
A fundamental problem is the input concerning the mobilities. Figure [9] demon- 
strates how the pulses of events starting at the outer mantle and observed in 
the core are stretched if a different set [19] of input parameters for the electron 
mobility is chosen. At Vu as = 3000 V and p imp = 0.66 ■ 10 10 /cm 3 , the pulse 
is 14% longer than with the default input parameters [18] - However, a simple 
time scaling makes the pulses indistinguishable. 

The length of a pulse is also strongly influenced by Pi mp . This is demonstrated 
in Fig.[10l At V Uas = 3000 V a change in p imp from 0.66 • 10 10 to 0.83 • 10 10 /cm 3 
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Fig. 9. Simulated core pulses for two different sets of input parameters for the 
electron mobility. The result of a time scaling of the longer pulse is also shown. 
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has an effect on the core pulse of similar size than changing the mobility. In 
this case, the shape of the pulse is also modified; a simple time scaling does 
not make the pulses inddistinguishable. However, the effect is quite subtle 
and will be difficult to observe in real data. Unfortunately, the value of pi mp is 
generally not known to very high precision. In addition, it often varies up to a 
factor of three along z. Therefore, it is unavoidable that pulses simulated for 
a particular detector or even detector section have to be scaled in time before 
more detailed comparisons to data can be made. Assuming a certain mobility, 
the length of the measured pulses can be used to derive the impurity density 
locally. 



9 Comparison to data 

The focus of the comparisons presented here is the drift of electrons. There- 
fore pulses as seen in the core electrode were used. Simulated pulses were 
compared to measured pulses produced by a well understood true-coaxial re- 
type detector [7J irradiated with a collimated Europium source. The 122 keV 
line was used to study energy deposits close to the detector outer mantle. In 
this case, the holes are quickly absorbed by the electrodes on this outer mantle 
and mainly the drift of the electrons is relevant for the shape of the resulting 
pulses. 

Two methods were used. Simulated pulses were either fitted to individual mea- 
sured pulses or to averages of sets of pulses obtained in a given configuration. 
The latter was done to reduce the effects of noise. 



9.1 Test environment 

Figure [TT] depicts the measurement setup. Three segments, numbered 13 to 
15, of the middle layer of the detector were scanned in at z = 0. A colli- 
mated 75 kBq 152 Eu source inside a copper collimator was used. The collimator 
pointed to the center of the detector and the spot had a 1 a diameter of about 
5 mm on the outer mantle of the detector. Details of the setup can be found 
elsewhere [7J. The operational voltage was 3kV. The step size of the scan 
was 5° in segment 14 and 10° in segments 13 and 15. The uncertainty in <p 
was approximately A0 = 2.5°. In total, 25 steps were made to cover 180° in 
azimuth. The pulses of the core and all segments were recorded. Events from 
the 122 keV 7-line were selected with a cut on the core energy of 2 a of the 
energy resolution, i.e. about ± 4keV. 

Simulations were done with the energy deposits located exactly on the outer 
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Fig. 11. Schematic of the setup used to scan three segments of an (60, 3z) segmented 
detector at z = 

mantle or distributed according to the interactions of the 122 keV photons 
from the 152 Eu source. A constant pi mp = 1.05 ■ 10 10 /cm 3 was implemented. 
This is the average of the values given by the manufacturer for the top and 
the bottom of the detector which differ by a factor of two. The electric field 
was calculated for a bias voltage of Vu as = 3 kV. The measured orientation of 
the crystal axes was used in the simulation; the (110) axis was set to 290°, 
almost aligned with one boundary of segment 15. 

An amplifier decay-time of 50 fis and a cut-off in bandwidth of 37.5 MHz 
were implemented according to the specifications of the electronics system. 
Electronic noise was not added to the simulated pulses to simplify direct com- 
parisons with individual measured pulses. 

9.2 Fits to individual pulses 

Simulated core pulses, C sim (t), were fitted to measured pulses, C meas , using 
three parameters only: the amplitude, A, the time offset, T , and the time 
scaling-f actor T scale : 

Cmeas {t) = A X C s i m (t /T sca i e + To) . (5) 

Figure Q~2] shows a randomly selected core pulse at = 200° and the fit of the 
pulse generated for this angle. This position falls between the (110) and the 
(100) axes. The actual values of the parameters A and T have no relevance; 
they are not simulated in real units. The important parameter is T sca \ e which 
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has a value of 0.90 ± 0.01. The xV dof = 184/146 = 1.3 for this fit is quite good. 
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Fig. 12. Fit of a simulated to a measured pulse at <j> = 200°. The dots represent the 
data with error bars indicating the noise level. 

Fig. [TBI shows the distribution of x 2 /dof values of all fits of pulses at <fi = 200°. 
The average % 2 /dof = 174/146 = 1.2. In general, the quality of the fits is very 
good. This indicates that the simulation describes the general shape of the 
pulses very well. To determine an overall T sca i e , a clean subset of measured 
pulses was selected by: 

• requiring x 2 < 220, i.e. x 2 /dof < 1.5, to eliminate background events which 
have probably more than one interaction in the segment; 

• rejecting events in which an intermittent DAQ problem might have affected 
the recording of the pulse. 

For each pulse, a value of T sca i e was determined. Figure [TJ] shows the distri- 
bution of the time scaling-factors of the pulses at <fi = 200°. The mean time 
scaling-factor of 0.91 ± 0.01 was extracted using a Gaussian fit. This was done 
for all positions in <fi. Figure [15] shows the (f> dependence of the mean time 
scaling-factors. The dotted vertical lines indicate the segment boundaries. A 
value of T sca i e = 1, independent of (p would be expected for a perfect simula- 
tion. However, considering the uncertainties in the input parameters, the level 
of agreement is quite astonishing. The time scaling-factors are systematically 
about 10% low. This means that the simulated pulses are too short, indi- 
cating that the velocities of the electrons are overestimated. As discussed in 
section [SI this could be due to too high mobilities assumed or due to a slightly 
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Fig. 13. x 2 distribution from the fits of the simulated to measured pulses taken at 
<fi = 200°. Only events with x 2 /dof < 1.5 are analysed further. The cut is indicated 
by a vertical line. 
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Fig. 14. Distribution of time scaling-factors from the fits of the simulated to te 
measured pulses taken at 4> = 200°. The result of a Gaussian fit is shown as a solid 
line. 



16 



underestimated pi mp . It is impossible to distinguish the two possibilities at 
this point. 



The (j) dependence was fitted with a straight line. The absence of an oscil- 
lation pattern with a period of 90° shows that the longitudinal anisotropies 
introduced in the model describe reality in a satisfactory manner. The obvious 
dependence of T sca i e on <p is difficult to explain. In germanium crystals, p imp 
often depends on r. If the central z-axis of the original crystal is not identical 
with the z-axis of the detector, a dependence of the observed average p imp 
can result. A further investigation whether this effect can explain the about 
5 % effect is under way. 
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Fig. 15. Mean time scaling-factors, points, as extracted for each <f>. The uncertainties 
were taken from the Gaussian fits. The result of a straight line fit to the points is 
given as a solid line. 



9. 3 Fits to averaged pulses 



It would be interesting to carefully study the detailed shapes of individual 
pulses to distinguish between different sources of deviation between simulation 
and data. However, Fig. IT2l shows that the noise level in this measurement does 
not allow this. In order to remove the effects of noise, the selected pulses at 
each scan position were averaged, even though the interaction points were 
spatially distributed according to the beam spot and the penetration depth of 
122keV photons. 
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Fig. 16. Fit of the averaged simulated to the averaged measured pulse at <j> = 200°. 
The statistical uncertainties in the averaged measured pulse are too small to be 
seen. 

The simulation was adjusted accordingly and the simulated pulses were also 
averaged. Averaged simulated pulses were fitted to averaged measured pulses. 
The result for <p = 200° is shown in Fig. US The fit yields T scale = 0.90 ± 0.01. 
This is compatible to the mean T sca i e = 0.91 ± 0.01 observed before. 

The x 2 /dof of the fit is ~ 8. This is quite large. In Fig. [T6J the deviations are 
also clearly visible. The simulated pulse has sharper edges most likely due to 
the overly optimistic cut of 37.5 MHz in the bandwidth; Figure [S] demonstrates 
how pulses get smoothed out by a cut in bandwidth of 10 MHz. At the very end 
of the pulse, the simulation shows an extra flattening. This would indicate that 
Pimp is slightly smaller than assumed. Nevertheless, the overall shape of the 
pulse is reproduced rather well and it will be tested soon, how well simulated 
pulses work as training sets for neural network based pulse shape analyses. 



10 Summary and Outlook 

A fully functional pulse shape simulation package for HPGe detectors was 
developed. It can easily be adapted to any common HPGe detector geometry. 
The results of the simulation are critically dependent on the input parameters. 
While some parameters are generic to germanium, the detailed properties of 
a detector like impurity concentration and axes orientation are also extremely 
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important. Simulations have to be adjusted to each detector. 



The longitudinal anisotropy of the elecron drift is correctly described by the 
simulation. Therefore the comparison of simulated to measured pulses can 
yield information about intrinsic detector properties like local variations of 
the impurity concentration. 

The quality of the simulated pulses will be further studied and the simulation 
will be used as input to pulse shape analyses. 
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